<html>
<!--
   Licensed to the Apache Software Foundation (ASF) under one or more
  contributor license agreements.  See the NOTICE file distributed with
  this work for additional information regarding copyright ownership.
  The ASF licenses this file to You under the Apache License, Version 2.0
  (the "License"); you may not use this file except in compliance with
  the License.  You may obtain a copy of the License at

       http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.
  -->
    <!-- $Revision: 924345 $ -->
<body>
<p>
This package provides classes to solve Ordinary Differential Equations problems
and also compute derivatives of the solution.
</p>

<p>
This package solves Initial Value Problems of the form
<code>y'=f(t,y,p)</code> with <code>t<sub>0</sub></code>,
<code>y(t<sub>0</sub>)=y<sub>0</sub></code> and
<code>dy(t<sub>0</sub>)/dp</sub></code> known. The provided
integrator computes estimates of <code>y(t)</code>,
<code>dy(t)/dy<sub>0</sub></code> and <code>dy(t)/dp</code>
from <code>t=t<sub>0</sub></code> to <code>t=t<sub>1</sub></code>,
where <code>y</code> is the state and <code>p</code> is a parameters
array.
</p>
<p>
The classes in this package mimic the behavior of classes and interfaces from the
<a href="../package-summary.html">org.apache.commons.math.ode</a>,
<a href="../events/package-summary.html">org.apache.commons.math.ode.events</a>
and <a href="../sampling/package-summary.html">org.apache.commons.math.ode.sampling</a>
packages, adding the jacobians <code>dy(t)/dy<sub>0</sub></code> and
<code>dy(t)/dp</code> to the methods signatures.
</p>

<p>
The classes and interfaces in this package mimic the behavior of the classes and
interfaces of the top level ode package, only adding parameters arrays for the jacobians.
The behavior of these classes is to create a compound state vector z containing both
the state y(t) and its derivatives dy(t)/dy<sub>0</sub> and dy(t<sub>0</sub>)/dp and
to set up an extended problem by adding the equations for the jacobians automatically.
These extended state and problems are then provided to a classical underlying integrator
chosen by user.
</p>

<p>
This behavior imply there will be a top level integrator knowing about state and jacobians
and a low level integrator knowing only about compound state (which may be big). If the user
wants to deal with the top level only, he will use the specialized step handler and event
handler classes registered at top level. He can also register classical step handlers and
event handlers, but in this case will see the big compound state. This state is guaranteed
to contain the original state in the first elements, followed by the jacobian with respect
to initial state (in row order), followed by the jacobian with respect to parameters (in
row order). If for example the original state dimension is 6 and there are 3 parameters,
the compound state will be a 60 elements array. The first 6 elements will be the original
state, the next 36 elements will be the jacobian with respect to initial state, and the
remaining 18 will be the jacobian with respect to parameters. Dealing with low level
step handlers and event handlers is cumbersome if one really needs the jacobians in these
methods, but it also prevents many data being copied back and forth between state and
jacobians on one side and compound state on the other side.
</p>

<p>
Here is a simple example of usage. We consider a two-dimensional problem where the
state vector y is the solution of the ordinary differential equations
<ul>
  <li>y'<sub>0</sub>(t) = &omega; &times; (c<sub>1</sub> - y<sub>1</sub>(t))</li>
  <li>y'<sub>1</sub>(t) = &omega; &times; (y<sub>0</sub>(t) - c<sub>0</sub>)</li>
</ul>
with some initial state y(t<sub>0</sub>) = (y<sub>0</sub>(t<sub>0</sub>),
y<sub>1</sub>(t<sub>0</sub>)).
</p>

<p>
The point trajectory depends on the initial state y(t<sub>0</sub>) and on the ODE
parameter &omega;. We want to compute both the final point position y(t<sub>end</sub>)
and the sensitivity of this point with respect to the initial state:
dy(t<sub>end</sub>)/dy(t<sub>0</sub>) which is a 2&times;2 matrix and its sensitivity
with respect to the parameter: dy(t<sub>end</sub>)/d&omega; which is a 2&times;1 matrix.
</p>

<p>
We consider first the simplest implementation: we define only the ODE and let
the classes compute the necessary jacobians by itself:
<code><pre>
public class BasicCircleODE implements ParameterizedODE {

    private double[] c;
    private double omega;

    public BasicCircleODE(double[] c, double omega) {
        this.c     = c;
        this.omega = omega;
    }

    public int getDimension() {
        return 2;
    }

    public void computeDerivatives(double t, double[] y, double[] yDot) {
        yDot[0] = omega * (c[1] - y[1]);
        yDot[1] = omega * (y[0] - c[0]);
    }

    public int getParametersDimension() {
        // we are only interested in the omega parameter
        return 1;
    }

    public void setParameter(int i, double value) {
        omega = value;
    }

}
</pre></code>
</p>

<p>
We compute the results we want as follows:
<code><pre>
    // low level integrator
    FirstOrderIntegrator lowIntegrator =
            new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10);

    // set up ODE
    double cx = 1.0;
    double cy = 1.0;
    double omega = 0.1;
    ParameterizedODE  ode = new BasicCircleODE(new double[] { cx, cy }, omega);

    // set up high level integrator, using finite differences step hY and hP to compute jacobians
    double[] hY = new double[] { 1.0e-5, 1.0e-5 };
    double[] hP = new double[] { 1.0e-5 };
    FirstOrderIntegratorWithJacobians integrator =
            new FirstOrderIntegratorWithJacobians(lowIntegrator, ode, hY, hP);

    // set up initial state and derivatives
    double     t0    = 0.0;
    double[]   y0    = new double[] { 0.0, cy };
    double[][] dy0dp = new double[2][1] = { { 0.0 }, { 0.0 } }; // y0 does not depend on omega

    // solve problem
    double     t     = Math.PI / (2 * omega);
    double[]   y     = new double[2];
    double[][] dydy0 = new double[2][2];
    double[][] dydp  = new double[2][1];
    integrator.integrate(t0, y0, dy0dp, t, y, dydy0, dydp);
</pre></code>
</p>

<p>
If in addition to getting the end state and its derivatives, we want to print the state
throughout integration process, we have to register a step handler. Inserting the following
before the call to integrate does the trick:
<code><pre>
    StpeHandlerWithJacobians stepHandler = new StpeHandlerWithJacobians() {
            public void reset() {}

            public boolean requiresDenseOutput() { return false; }

            public void handleStep(StepInterpolatorWithJacobians interpolator, boolean isLast)
                throws DerivativeException {
                double   t = interpolator.getCurrentTime();
                double[] y = interpolator.getInterpolatedY();
                System.out.println(t + " " + y[0] + " " + y[1]);
            }
    };
    integrator.addStepHandler(stepHandler);
</pre></code>
</p>

<p>
The implementation above relies on finite differences with small step sizes to compute the
internal jacobians. Since the ODE is really simple here, a better way is to compute them
exactly. So instead of implementing ParameterizedODE, we implement the ODEWithJacobians
interface as follows (i.e. we replace the setParameter method by a computeJacobians method):
<code><pre>
public class EnhancedCircleODE implements ODEWithJacobians {

    private double[] c;
    private double omega;

    public EnhancedCircleODE(double[] c, double omega) {
        this.c     = c;
        this.omega = omega;
    }

    public int getDimension() {
        return 2;
    }

    public void computeDerivatives(double t, double[] y, double[] yDot) {
        yDot[0] = omega * (c[1] - y[1]);
        yDot[1] = omega * (y[0] - c[0]);
    }

    public int getParametersDimension() {
        // we are only interested in the omega parameter
        return 1;
    }

    public void computeJacobians(double t, double[] y, double[] yDot, double[][] dFdY, double[][] dFdP) {

        dFdY[0][0] = 0;
        dFdY[0][1] = -omega;
        dFdY[1][0] = omega;
        dFdY[1][1] = 0;

        dFdP[0][0] = 0;
        dFdP[0][1] = omega;
        dFdP[0][2] = c[1] - y[1];
        dFdP[1][0] = -omega;
        dFdP[1][1] = 0;
        dFdP[1][2] = y[0] - c[0];

    }

}
</pre></code>
With this implementation, the hY and hP arrays are not needed anymore:
<code><pre>
    ODEWithJacobians ode = new EnhancedCircleODE(new double[] { cx, cy }, omega);
    FirstOrderIntegratorWithJacobians integrator =
            new FirstOrderIntegratorWithJacobians(lowIntegrator, ode);
</pre></code>
</p>
</body>
</html>
